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Abstract 

Epitaxial self-assembled quantum dots (SAQDs) are of interest for nanostructured optoelectronic and electronic 
devices such as lasers, photodetectors and nanoscale logic. Spatial order and size order of SAQDs are important 
to the development of usable devices. It is likely that these two types of order are strongly linked; thus, a study of 
spatial order will also have strong implications for size order. Here a study of spatial order is undertaken using a linear 
analysis of a commonly used model of SAQD formation based on surface diffusion. Analytic formulas for film-height 
correlation functions are found that characterize quantum dot spatial order and corresponding correlation lengths that 
quantify order. Initial atomic-scale random fluctuations result in relatively small correlation lengths (about two dots) 
when the effect of a wetting potential is negligible; however, the correlation lengths diverge when SAQDs are allowed 
to form at a near-critical film height. The present work reinforces previous findings about anisotropy and SAQD order 
and presents as explicit and transparent mechanism for ordering with corresponding analytic equations. In addition, 
SAQD formation is by its nature a stochastic process, and various mathematical aspects regarding statistical analysis 
of SAQD formation and order are presented. 

1 Introduction 

Epitaxial self-assembled quantum dots (SAQDs) represent an important step in the advancement of semiconductor 
fabrication at the nanoscale that will allow breakthroughs in optoelectronics and electronics. fl]|2]E|4][5]|6l|7l|8][9] 
ITOl [TT1 fl2ll Most frequent optoelectronic applications are high efficiency lasers with exotic wavelengths or photode- 
tectors. d (3] El E] IS] IZl El [TT] [121 SAQDs are the result of a transition from 2D growth to 3D growth in strained 
epitaxial films such as Si x Ge 1 _ 2 ,/Si and Inj-Ga^jjAs/GaAs. This process is known as Stranski-Krastanow growth or 
Volmer- Webber growth. fl3l IT1 Il4l [El . 

In applications, order is a key factor. There are two types of order, spatial and size. Spatial order refers to 
the regularity of SAQD dot placement, and it is necessary for nano-circuitry applications. Size order refers to the 
uniformity of SAQD size which determines the voltage and/or energy level quantization of SAQDs. It is reasonable 
to expect that these type of order are linked, and it is important to understand the factors that determine SAQD order. 
Further understanding should help in the design and simulation of both spontaneous "bottom up" self-assembly and 
directed or guided self-assembly to enhance SAQD order. Ifl6l [T71 [181 [T9l l20l |2T1 122[ [23]Here, an elaboration of and 
further application of a linear analysis of SAQD order [24 J is presented. The work reported here forms the basis of a 
non-linear theory and modeling of SAQD order that will be reported in future work. 

In [24 1 it was reported that one could calculate a correlation function using a linearized model of SAQD formation. 
This correlation function included two correlation lengths that could be used to describe SAQD order. It was also found 
that one effect of a hypothesized wetting potential was to enhance SAQD order when growth occurs near the critical 
film height for 3D growth. Here, these results are expanded to create a more rigorous linearized theory of SAQD order 
that will inform non-linear theories. In particular, the model is generalized to any model that combines local energy 
effects such as surface energy density and non-local elastic destabilization, and the procedure for predicting order 
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based on any linear theory with peak wavelengths is presented. The hypothesized effect of elastic anisotropy in [24] 
is verified with calculations using linear anisotropic elasticity theory. l25l [26 1 Details such as statistical fluctuation 
and convergence are also addressed along with a discussion of the possible forms of linear anisotropic terms in SAQD 
growth kinetics, and the effect of an atomic-scale cutoff in the continuum theory is addressed. Finally, the order 
enhancing effect of growing near the critical threshold is explored in more detail using calculations appropriate to 
Ge/Si SAQDs. 

In the literature, two modes of SAQD formation are generally discussed, the thermal nucleation mode and the nu- 
cleationless mode. Il27l l28l l29ll In the thermal nucleation mode, a 2D film surface is metastable, and the formation of 
individual quantum dots is thermally activated. l27l . This growth mode leads to the formation of individual quantum 
dots as uncorrelated or loosely correlated discrete events at essentially random locations. In the nucleationless mode, 
the 2D film surface transitions from stable (or metastable) to unstable. In this mode, dots form everywhere at once 
appearing at first as a cross-hatched ripple-like disturbance on the 2D film surface and then maturing into recogniz- 
able individual dots. [27, 30, 28, 3T1 l32ll These two modes are probably connected via an encompassing conceptual 
and mathematical modeQ and perhaps some of what is observed experimentally is in fact a hybrid mechanism. In 
agreement with intuition, it appears that the nucleationless mode leads to a more ordered dot pattern than the thermal 
nucleation mode that is dominated by randomness.]^] Thus, the presented analysis applies to the nucleationless mode. 

There are various implementations of nucleationless growth models [28 37, 38, 39, 40, 18, 34|, although, there is 
also a great deal of commonality among these models. In particular, they all include a non-local elastic effect and local 
surface energies and/or local wetting energies. Here, a linear analysis of quantum dot order resulting from this class of 
model is presented. Particular note is taken of the effects of stochastic initial conditions crystal anisotropy in general, 
elastic anisotropy in particular, and the effect of varying film height as a control parameter as first introduced in ll33l . 
A simple model similar to l28l l37l [381 l40l [T8l is presented to produce numerical examples and explore the effects 
of the average film height. Concurrently, a more abstract and general model is presented and analyzed that includes 
non-local elastic strain effects, and a local combined surface and wetting energy. The linear model with stochastic 
initial conditions and deterministic film height evolution will pave the way for more sophisticated analysis involving a 
non-linear model of stochastic film height evolution. 

As previously stated, one of the goals in the present work is to further explore the role of the wetting poten- 
tial during growth near the stability threshold in film height. A wetting potential has been included in the analysis 
and simulations in [38, 33] EJUS). Although somewhat controversial, the wetting potential plays an important phe- 
nomenological role. It ensures that growth takes place in the Stranski-Krastanow mode: that a 3D unstable growth 
occurs only after a critical layer thickness is achieved, and that a residual wetting layer persists. The physical origins 
and consequences of the wetting potential are discussed in flTll28l . The analysis presented here is usable in models 
that neglect the wetting potential by simply setting it to zero. Another possibility is simply that the wetting potential 
is simply an approximation to the stabilizing effect of intermixing. [42] That said, if the wetting potential is real, the 
present analysis shows that it is beneficial to SAQD order to grow near the critical layer thickness. 

The presented analytic formulas and linear analysis are intended to complement existing numerical models of 
SAQD order. l43l [37l l44l l45l and to form a basis for future non-linear analytic analysis of SAQD order. The current 
findings agree with previous work on the beneficial effects of elastic anisotropy to enhance in-plane order. 

The linear analysis, of course, represents a simplification of the film evolution, and it applies only to the initial 
stages of SAQD formation when the nominally flat film surface becomes unstable and transitions to three-dimensional 
growth. However, the small surface fluctuation stage of SAQD growth determines the initial seeds of order or disorder 
in an SAQD array; thus, the small fluctuation stage should have an important influence on the final outcome. At later 
stages when surface fluctuations are large, there is a natural tendency of SAQDS to either order or ripen [33, 37ll46l 
[39]|47| Ordering systems tend to evolve slowly due to critical slowing down |39], while ripening tends to diminish 
order further. fl37l Thus, it is possible that the linear model could, in fact, yield good predictions of SAQD order. 
The simplification and linearizion facilitates the development of analytic solutions that are most transparent, easily 
portable to multiple material systems and have no effective limit on system size. Finally, it is virtually impossible to 
have a thorough understanding of the full non-linear model without first having a thorough understanding of the linear 
behavior. 

The remainder of the paper is organized as follows. Section[2]presents the physical assumptions and mathematical 

'it is likely that there a transition from stable, to metastable and finally to unstable. The analysis presented in |33| would appear to support such 
a view where the film height acts as the control parameter driving the transition. There is also some controversy regarding whether all dot growth is 
nucleationless or not. 13411321 

2 Compare various figures in I29II35IIT41 [3111361 . 
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approximations used to model film growth. Section [3] discusses the stochastic initial conditions and the resulting 
correlation functions and correlation lengths. SectionPfl presents a procedure for estimating SAQD order with an 
application to Ge dots on a Si substrate. Section [5] presents conclusions, while Appendices |AfF| present additional 
calculational details. 



2 Modeling 

The formation of SAQDs is modeled as a deterministic surface diffusion process with stochastic initial conditions. The 
resulting equations and ultimately the sought after correlation functions are different depending on whether the film 
surface is treated as one-dimensional isotropic, two-dimensional isotropic or two-dimensional anisotropic. The ID 
and 2D isotropic cases are discussed first, and then the essential differences of the 2D anisotropic model are presented. 
The stochastic initial conditions need to be expressed in terms of the correlation functions that are also use to analyze 
order; consequently, the discussion of the initial conditions is deferred to Sec. |3.2| 

It should be noted that the results presented here are fairly general. There has been a good deal of recent work 
refining the modeling of nucleationless growth processes to incorporate various phenomenological aspects of SAQD 
growth. For example, the inclusion of orientation-dependent surface energy fl38l . strain-dependent surface energy [34] 
and explicit modeling of atomic species segregation and film-substrate inter diffusion. Il48l Two models are presented 
here. One is a simple concrete example. It is the simplest model one can use including elastic effects surface energy 
and wetting energy. The second model is more abstract and describes the general case of a local potential energy 
that depends on both the film height and film height gradient. One effect that is not examined here is that of mixed 
4-fold and two-fold symmetry. Such a mixing can occur due to diffusional anisotropy or surface energy anisotropy. 
(Sec. |2.2.1.2| and Appendix [D]). However, a similar analysis procedure should work for these cases as well. The general 
procedure for possible application to other models is discussed in Sec. |3.5| 

The following discussion will use abstract vector notation, e.g. k instead of ki, etc. Also, because it is sometimes 
computational expedient to perform one -dimensional modeling [ 24 39j [T7] |42j] •. the case of a one dimensional surface 
with two dimensional volume is discussed along with the case of an isotropic 2D surface. To facilitate this combined 



discussion, the dimensionality of the surface will be denoted as d. In Sees. 3.3 and 3.4 d — 1, 2 will be substituted as 
appropriate. Finally, much of the calculation involves reciprocal space. The convention used for the Fourier transforms 
is 

/(x) = J Ae !k '7 k , and f k = (2n)- d J d d x e -* kx /(x) 

following the example of [ 



2.1 ID and 2D Isotropic model 

This discussion pertains to both the ID model and the 2D isotropic model. The formation of SAQDs is modeled as a 
surface diffusion process where the film height is a function of the lateral position and time. The system is treated as 
deterministic with stochastic initial conditions. First, the general non-linear governing equations are presented. Then, 
the linearized form is presented. Finally, the key behavior is reviewed. 

The mathematical model uses film height, 7i(x,t) as the dependent variable and the horizontal position x and 
time t as the independent variables. The film height evolves over time due to surface diffusion driven by a diffusion 
potential /i(x, t) and a flux of new material Q. The surface velocity is thus 

v n =n z d t n=-Vs-VV s ^,t) + Q (1) 

where n z is the vertical component of the surface normal n z = [1+ (V7i) 2 ] -1 / 2 , Vs is the surface gradient, T> is the 
diffusivity, and Vs- is the surface divergence. 



2.1.1 Energetics 

The diffusion potential /i(x, t) must produce Stranski-Krastanow growth. Thus, it must contain an elastic term that 
destabilizes film growth, a surface energy term that stabilizes planar growth and a wetting energy that ensures a wetting 
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layer. The diffusion potential can be derived from a total free energy. 



J~ J~ slast ~t~ surf. 

dVuo- 

volume 



surface 



eL4 surf .7 + / dAW(H) 



where uj is the elastic energy density, 7 is the surface energy density, W(Ti) is the wetting energy density. The 
last integral corresponds to JF we t, and whether the integral should be taken over the film surface or the substrate is 



ambiguous. The "simple" model (Sec. 2.1.1.1 1 assumes that the integral is over the substrate, while the "general" 
model (Sec. |2. 1 . 1 .2| > can accommodate both cases. 

2.1.1.1 simple form The simplest possible model results if the integral corresponding to T we t is taken over the 
lateral positions x rather than over the actual free-surface. In concrete terms, one can use dV = d 2 x.dz and d^4 S urf. = 



l + (VW(x))" 
T = 



1/2 



to obtain the expression 

<i 2 x<izti>[7i](x, z) - 



volume 



dh 



x-plane 



1 + (VW(x))* V2 7+ W(H(x)) 



(2) 



where the indicates that the elastic energy density is a non-local functional of the film height, H. The diffusion 

potential ji can be found, similar to IPT51 . by differentiating T with respect to the surface motion (Appendix |A.1|>, 



/x(x) = Q8!F/S?i(x). Doing so for Eq. Q (Appendix A. 2 1, 

fi{x) = n [u(x) - 7«(x) + W (W(x))] 



(3) 



where Q, is the atomic volume, cj(x) is the elastic energy density at the film surface (implicitly uj[H] (x,7i(x))), 

is the total surface curvature, and W (H) — dn(x) W (W(x)) is the 



k = V- |VW(x) [l + (VH(x)) 
derivative of W (7i(x)) evaluated at x. 

2.1.1.2 general form It should be noted that Eq. ([3} is not the same diffusion potential used in [38]. The wetting 
potential used there can be derived by taking W(H) as an energy density of the free surface, not a density in the 
x-plane. Expressions like Eq. ([3]) and Eq. (1) in [38 1 are part of a larger class of surface evolution models with more 
or less the same linear behavior. 

The surface and wetting energy can be combined and incorporated into a more general form, with a total free 
energy T sw and a free energy density F SW (H., VW) that depends on the film height W(x) and the film height slope or 
orientation VH(x), The total free energy is thus 



-l2 



(4) 



d -Kdzuj[H] (x, z) 



volume 



d 2 *F sw (H(x),VH(x)) 



x— plane 



F sw may not necessarily be the sum of separate surface energy and wetting energy contributions. It need only be a 
local function of H and VW. The corresponding diffusion potential is 



Kx) = n L(x) + ^°)(x) - v • f(°J)(x) 



(5) 



where Fsw^ indicates the m th derivative with respect to TL and the n th derivative with respect to V7i . Fgffl (x) = 
9h(x)Fsw (H(x),VH(x)) and each vector component of F^(x) is F^^(x) = 9[vw(x)].-^»io ("^( x ) 5 V^( x )) 
One can obtain the results of the simple model (Eqs. |2]i and <|3j) by setting 



p — 



1/2 

l + (VH(x)) z 7 + W(W(x)) 



(6) 



A diffusion potential like Eq. (1) in [38 1 can be obtained by setting 



,2 1 ' 2 



(VW(x)) z [ 7 (VW(x)) + W(W(x))] 
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This is different from Eq. (|6) in two ways. First, the surface energy density depends on the surface orientation. Second, 

r 21 1/2 

the Jacobian, J = 1 + (V7i(x)) multiplies both the surface energy density and the wetting potential. Despite 
these differences, the common form of the diffusion potential (Eq. (|5]l) among different models suggests that they 
might all lead to similar linearized forms and behavior. 

2.1.1.3 Linearization The diffusion potential is now linearized about the average film height TL. In general, one 
can control the amount of deposited material, and thus the average film height TL. It is therefore useful to decompose 
7i(x) into the spatially averaged mean value and fluctuations about the average. Similar to 11281 . 

H = H + h(x,t). (7) 

In the present calculation, TL is specified as constant in time. This assumption corresponds physically to a fast deposi- 
tion and then an anneal. It is not too difficult to generalize to a time dependent TL, but that is beyond the scope of this 
manuscript. In 1 38 . 49 1, deposition and evaporation is explicitly modeled. 

All terms in /x(x, t) are now kept to only linear order in h(x,t). The elastic energy density oj is a non-local 
functional of ft,(x, t) |40|; however, the equations generating w(x) are translationally invariant. Thus, it is convenient 
to use reciprocal space for the linearization. The curvature is trivially linearized as k(x) — > V 2 /i(x) in real space or 
Kk — ► — fc 2 /ik in reciprocal space. The linearized elastic strain energy ui can be found in reciprocal space as in ifTSIl 
to be Wk = — 2M(1 + i/)e 2 n /ik, where M = E/(l — v) is the biaxial modulus, E is the Young modulus, v is the 
Poisson ratio, and e m is the film-substrate mismatch strain. This formula neglects possible differences in elastic moduli 
between the film and substrate as in [28 1, but a similar method of analysis should apply to that case as well. Linearizing 
Eqs. ([3} and ([5} in reciprocal space, /ik is proportional to with a proportionality coefficient that depends on k and 
TL. 

Miin,k = /(k,W)/i k (8) 

where /(k, TL) for three different isotropic cases, corresponding to Eqs. ^ and ([5]), and an abstracted general form, is 
given by 

( fl [-2M(1 + v)t 2 m k + 7 fc 2 + W"{TL)} ; case a (Eq. ^) 
/(k, TL) = \ fl [-2M(1 + v)e' 2 m k + F 02 fc 2 + F 20 ] ; case b (Eq. ^) . (9) 
I— ak + bk 2 + c ; case c (general) 

Due to isotropy, /(k,7Y) is independent of the direction of k, and only the wave number, k — ||k||, appears in the 
right hand side. is the second derivative of F sw with respect to TL, and Fsffl the second derivative of F sw with 

respect to V7i. Fg^ and Fgffl depend on TL only; thus they are constants in the present analysis. See Appendix 
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for more precise definitions and the derivation of /(k, TL). Using Eq. |6|, produces Fgffl = 7 and Fg^ = W"(TL) 
which is identical to the simple case of Eq. a. Case c, labeled as "general" where a, b, and c depend implicitly on 
TL shows that /(k, TL) for cases a and b have the same relatively simple form. It also emphasizes the dynamic effects 
as opposed to the physical causes. There is a destabilizing term, — ak, a short wavelength cutoff term, bk 2 , and a term 
that stabilizes the entire spectrum, c. 

Despite the label "general," there are of course limitations to the application of Eqs. ([8| and Q. For example, 
there has been recent work on the effects of strain-dependent surface energies. [34] The second form can not represent 
such an effect because the derivation assumes that the surface energy only depends on local quantities, {TL and V7i) 
whereas the strain effect is non-local. However, it is reasonable to conjecture that a more detailed analysis of the 
effects of a strain dependent surface energy term would produce a coefficient function /(k, TL) not very different from 
the case c "general" form of Eq. (|9). Thus, the following analysis may very well apply to this more exotic model, but 
more study is needed to be certain. 



2.1.2 Dynamics 



As discussed in Sec. 2.1.1 the dynamics are derived assuming no flux of new material (Q = 0) and keeping only 
terms to linear order in the height fluctuation, /i(x,t). Under these assumptions, Eq. ([T| can be decomposed into a 
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Table 1 : Characteristic wave-numbers, characteristic times and associated dimensionless variables for the three cases 
addressed in Eq. |9). 





k c 


tc 


a 


(3 


case a 




7 3 


k/fc c 


jW"(H) 
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16'DS2M 4 (l + ^) 4 e8 1 


4M 2 (l+i/) 2 £ 4 n 


case b 


2M{l+v)e 2 m 






p (02) p(20) 


p(02) 




4M J (l+i/) ! eJ 1 


case c 


a/b 


6 3 /(Ma 4 ) 


k/k c 


cb/a 2 



trivial equation for H and an equation for the film height fluctuation by inserting Eq. |7}. 

dU/dt = (10) 
dthfr) = -V ■ PV Win (x) (11) 

where /iu n (x) is the inverse Fourier transform of Eqs. <|8j and Q, and it depends implicitly on the average film height 
H. Note that the time dependence is implicitly while the coordinate dependence is explicit. The explicit coordinate 
dependence serves to distinguish Assuming that the diffusivity T> is constant, the Fourier transform of Eq. ( fTT| gives 
the linearized differential equation for the evolution of each Fourier component. 

d t h k = -Vk 2 ^ = -Vk 2 f(k,H)h k . (12) 

Solving Eq. ([12), 

&k(t) = M0)e CTkt ; (13) 

<r k - -Vk 2 f(k,H). (14) 

The surface evolves in reciprocal space as an initial condition, /ik(0) multiplied by an envelope function, e CTk *. For 
most values of Ti, this envelope function has a peak. As time passes, this peak narrows and can be approximated by a 
gaussian. To analyze this behavior, appropriate dimensionless variables are defined. Then, the stability of the film is 
discussed. Finally, is expanded about its peak to aid analytic calculations. 

The time dependent behavior of the film height fluctuations is facilitated by using a characteristic wave number, 
characteristic time and related dimensionless variables. For the "general" case c of Eq. (|9), the characteristic wavenum- 
ber is k c — a/b, and the characteristic time is t c = 1/ (T>Q,bk 4 ) = 6 3 / (Vila 4 ) . These characteristic dimensions can be 
used to define a dimensionless wave vector, a. = k/fc c and a dimensionless wetting parameter j3 = c/ (bk 2 ) = cb/a 2 . 
One can also define a dimensionless time, r = t/t c . To obtain the corresponding characteristic scales for cases a and 
b, one merely has to plug in the appropriate substitutes for a, b and c and follow the pattern. For example, for case 
a, make the substitution a — * H2M(2 + v)^, etc. Table [T| summarizes these values for all three cases. For all three 
cases, /(k, Ti) and the growth constant er^ reduce to the following forms: 

/(k,W) = f(k c a,n) = nbk 2 c (-a + a 2 + l3) (15) 
^k = Ok c _ a = t~ x u 2 (a -a 2 - 0) , (16) 

where a = \\a\\ = k/k c is the dimensionless wave number. These forms are plotted in Figs.[lJ and[2] Fig.[jJ shows 
f{\s.,1i.)/Vtbk 2 vs. a for an isotropic or one dimensional surface. Figs. [2] shows t c Ok vs. a for a 2D anisotropic 
surface (Sec. |2.2) . However, the curves marked 0° are identical to the dispersion relation for a ID or 2D isotropic 
surface (compare Eqs. (|9) and (|23)). 

2.1.3 Peaks 

The peak growth rate and the corresponding wavenumber k can be found from Eq. ([16). has a peak at ko — k c ao 
where 

a = I (3+^9-32/3) . (17) 
Expanding <7k about this peak to second order in k — ko, 

fk ~ c - \ a 2(k - k ) 2 
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Figure 1: Dimensionless diffusion potential pref actors vs. dimensionless wave number, (a) The one dimensional or 
isotropic case with (3 = 0.3. (b) The elastically isotropic case with anisotropy = 0.1 (see Eq. (|22]>). 
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(a) Growth Constant vs. Wave Number 
0.015 F 
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(b) Growth Constant vs. Wave Number 

0.1 




Figure 2: Dimensionless growth constant vs. dimensionless wave number. Curves are plotted for the elastically 
anisotropic case, but the curves marked 0° are the same as for the isotropic cases. In (a), /3 = 0. In (b) [3 = 0.2. 
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Figure 3: Exponential Envelope e CTk ' as function of a for j3 = 0.208 and t/t c = 100. (a) 2D isotropic surface, (b) 2D 
anisotropic surface with e = 0.1236. 

The two constants are 

o~o = -C^o («o - 2/3) , (18) 

and 

0-2 =t~ 1 k~ 2 (3a -4p). (19) 
Inserting this approximation for into Eq. ( fT3} , 

fc k (t) = /i k (0)e CTot e~5' T2t( ' £ " fco)2 . (20) 

The individual initial surface fluctuation components grow with a gaussian shaped envelope. An example of this 
envelope is plotted in Fig. [3|a). Notice that in two dimensions, the envelope forms a ring as the peak is about the 
wave-number fco but not about any particular point in the k-plane. 

2.1.4 Stability and wetting potential 

Stranski-Krastanow growth is marked by a transition from stable two-dimensional growth to unstable three-dimensional 
growth once a critical height 7i c is reached. [1] Eqs. ( (IT) , ([18) and ( |20| i are useful for analyzing the transition from 
stable to unstable growth. In order for this transition to occur, there must be some stabilizing term in the diffusion 
potential. In the present model, this means that there must be some surface energy-like term that varies strongly with 
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film height. This condition equates to stating that W"(n) or F^® or c (Eq. must be rather large if Ti < Tie- 
However, as H increases, these terms are reduced. Finally, when H > H c , this term is no longer capable of stabilizing 
the film against fluctuations of all possible wavelengths. 

The critical value TL C can be found using the analysis from [33 1. By inspection of Eqs. ([SJ, (|9]l and ( [12) , modes 
with / > increase the total free energy T as they grow; thus, they are stable and decay with time. Modes with 
/ < decrease the total free energy T as they grow; thus, they are unstable and grow with time. This growth and 
decay rule is easily verified by inspection of Eq. ( |14) , Thus, stable growth occurs when /(k, TL) > for all values of 
k, and unstable growth occurs when /(k, H) < for some values of k. Thus, the transition from stable to unstable 
growth occurs when the minimum value of /(k, H) just becomes negative. Using the same dimensional analysis as in 
the previous section and following the discussion of 1551 . one finds that the minimum value, / m ; n = ilbk 2 ((3 — 1/4), 
occurs at k m i n /k c = a m ; n = 1/2. / m ; n first becomes negative, and the transition to unstable growth occurs when 
the dimensionless wetting parameter (Table [TJ drops to a critical value, (3 = 1/4 . /3 > 1/4 stable 2D growth, and 
(3 < 1/4 unstable 3D growth. It is reasonable to suppose that W^H), W"(H), and thus are positive monotonically 
decreasing functions of H so that the interface becomes less important for large values of H. For example, in [50| it is 
assumed that W(H) = B/H, where B is constant. When (3 — > 0, corresponding to large H, the case discussed in [28] 
is obtained. A similar analysis can be done for cases b and c once one specifies how the terms Fsffl and Fgffi or a, b 
and c depend on H. 

Using a guessed form for a wetting potential, one can find the critical film height H c by setting (3 — 1 /4 . Applying 
this condition to case a in Eq. ([5J 

W"(H C ) = jk 2 c /A. 

Using the wetting potential of ll50l as an example, W(TL) — B/H, 

H c = {/8B/{^kf) = </82?7/(2M(l + I /)e2 i )2. (21) 

Conversely, one can fit a wetting potential to an observed or reasonable critical layer thickness from the same condition. 
Using the example wetting potential from [50], 

{2M{l + v y m fHl 
B - 8^ ' 

as stated in (5010 



2.2 2D Anisotropic case 

Crystal anisotropy leads to a dispersion relation that is both quantitatively and qualitatively different from the 
isotropic case. Here the effect of elastic anisotropy is discussed in most detail. Other sources of anisotropy are 
the surface and wetting energies. For example, in l!58l the surface energy density is orientation dependent which 
introduces a possible anisotropy in the dispersion relation. Possible sources of anisotropy are an anisotropic elastic 
stiffness tensor, an orientation dependent surface energy or wetting potential or anisotropic diffusion. As discussed 
below, the form of anisotropy to linear order in the height fluctuation, h, is somewhat restricted. Results are presented 
for 4-fold symmetric surfaces, that is surfaces that have invariant dynamic evolution laws when rotated by 90° . Possible 
complications arising from 2-fold symmetric anisotropic terms (with 180° rotational symmetry) are also discussed. As 
for the isotropic case, first the energetics are discussed, then the dynamics, and finally the expansion about the peaks 
in the dispersion relation, ov 

2.2.1 Energetics 

The discussion of energetics will first treat the effects of elastic anisotropy and then anisotropy resulting from surface 
or wetting like terms. 

3 This result from 1501 corresponds to the choice F,n„ CH. VW = \l 1 + ( VH) 2 7 + W{'H). However, the numerical model in 1501 appears to 

useF sw (H, VW) = yl + (VW) 2 b(VW) + W{H)}. This difference should lead to a slightly different critical film height in their numerical 
model from the one that they predicted (Eq. \2\\ ). 
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Elastic Energy vs. Direction 




'■ 7 o 



10° 20° 30° 

6. (degrees) 



Figure 4: Plot of£g k /(Me^ n )for various materials. Symbols indicate values calculated using Appendix|c] Solid lines 
are the interpolation (Eq. ( p2| ) using the values from Table [2] 



Table 2: Elastic c onstants [51 1 and calculated values (see Appendix [C| for various materials of inte rest at T = 300K. 





Cll 

io 11 ^! 


Cl2 

10 n ^| 

Ctrl' 5 


c 44 
lO 11 ^ 

cm- 5 


M 
10 n ^| 


£ ° 

Mel 


£450 

Mel 




Ge 


12.60 


4.40 


6.77 


13.93 


2.16 


1.906 


0.1176 


Si 


16.60 


6.40 


7.96 


18.07 


2.22 


1.997 


0.1005 


InAs 


8.34 


4.54 


3.95 


7.94 


2.70 


2.09 


0.226 


GaAs 


11.90 


5.34 


5.96 


12.45 


2.15 


1.87 


0.1302 



2.2.1.1 Elastic anisotropy One would like to obtain a simple symbolic expression for the elastic energy density at 
the free surface, w^, to first order in h* for the elastically anisotropic case. Similar discussions can be found in Il25ll26l . 
For the isotropic case, cu* = — 2M(1 + z/)e^ l /ii c . For the anisotropic case, 

Wk = -£e k kh k 

where the prefactor £g k is the decrease in elastic energy at the surface per unit wave number (k — ► 1) and unit amplitude 
(/ik — > 1) . It is not constant, but instead depends on the (?k, the angle that k makes with the a;— direction. The case of a 
cube-symmetry elastic stiffness tensor such as for Si is considered where one must specify three elastic constants en, 
C12 and C44. iBTI . Growth on a (100) surface will produce an elastic energy prefactor £g k that is four-fold symmetric 
(symmetric upon rotations by 90°). A procedure similar to Il25l 1251 based on a first order perturbation analysis is 
followed (Appendix O. A relatively simple interpolation formula 11241 is hypothesized and then verified numerically. 

The interpolation procedure, suggested in ll24l uses the lowest possible order expansion in sin(#k) and cos(#k) 
that has the appropriate four- fold symmetry and then interpolates between 6*k = 0° and 6>k = 45°. Thus, 



£ 9k =£ o (l-e A sm 2 (20*)) 



(22) 



where ea = (£ ° — £45°)/£o° is an anisotropy factor. This lowest order form turns out to be a very good fit to 
numerical calculations (Fig. [4]). Table [2] gives values of £ ° and €a for some systems of interest. In the elastically 
isotropic case, £0° — £450 = 2M(1 + v) so that ea = 0. 

There are two important differences from the elastically isotropic case. The first is obvious, that £g k depends on 
angular orientation, 9*. The second is that the peak value of ui* is not the same as that for the elastically isotropic case 
because in general, £q° ^ 2M(1 + v). In lF24l . where the purpose was simply to investigate the mechanism by which 
elastic anisotropy effects order, this second difference was neglected. 
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2.2.1.2 Surface and Wetting Energy Anisotropy The surface energy and wetting potential can be additional 
sources of anisotropy if they depend on the surface orientation so that 7 — * 7(V7i) or W(Ti) — > W(H, V7i) (for 
example, 115211381 ). Then, to first order in h , 

Msurf.,k = (jk 2 + k • 7" • k) /i k 

where 7" is the (2 x 2) matrix or Hessian matrix that results from taking the second derivatives of r y( V7i) with respect 



to the two components of V7i (Appendix B.l 1. Similarly 



wet 



k = ft(V( 20 > + k.w(° 2 > -k) /i k 



For 



where VK (20) and W (02) are the second derivatives of W(H, VH) with respect to H and VH (Appendix B.l 
both /isurf^k an d Mwet.k. the first term is isotropic, and the second term contains any possible anisotropy. 

The rank of the 7" and W( 02 ) matrices greatly restricts the possible forms of the additional anisotropy. These 
(2 x 2) matrices must be either two-fold symmetric or perfectly isotropic. Thus, if the surface energy and wetting 
potential are four-fold symmetric as £g k is, then 7" — > 7", a scalar, and W"( 02 ) — > W^ 02 \ a scalar, and neither one 
contributes any additional anisotropy. They do, however, help to stabilize or further destabilize the 2D surface as they 
add terms proportional to k 2 . The effect of these additional terms is indistinguishable from the effect of varying the 
value of the surface energy density, 7. 152, 3D 

It should be noted that the (100) surface of a diamond or zinc-blend structures allows for anisotropy that is only 2- 
fold symmetric (rotations by 180°). Thus, they could "break" the four-fold symmetry that occurs when one considers 
the elastic anisotropy alone. However, this "broken" symmetry is somewhat dubious because even the diamond and 
zinc -blend structures have a screw symmetry (rotations by 90° and translation in the [100] direction by half a lattice 
vector). Thus, if for example, W(H, V7Y) is anisotropic with two-fold symmetry to linear order, there must be a 
fast oscillation with changes in the film height Ti. In Appendix [Dj a similar term related to anisotropic diffusion is 
discussed. There does not appear to be any evidence for this two-fold symmetry in the case of (100) surfaces of IV 71 V 
systems such as Ge/Si, but in III-V/III-V systems the four-fold symmetry of the (100) surface may indeed be "broken" 
in this way corresponding to either a surface energy anisotropy or a diffusional anisotropy. P33ll54l . Further analysis 
of such terms in any more detail would greatly complicate the present discussion, so it is left for future work. Most 
of the modeling literature avoids this complication by not including the symmetry-breaking of the zinc-blend surface, 
for example (25] |26] [38]]. 

One can perform a similar analysis of the combined surface and wetting potential, F sw (TC, V7i ) (case b). To linear 



order the resulting anisotropic diffusion potential is (Appendix B.2 1 



Again, Fiu? is a rank 2 tensor, and all of the same symmetry considerations apply here as well. 

Because the two-fold symmetry anisotropic terms are excluded from the current discussion, and isotropic terms 
simply "renormalize" the effective of surface energy, there will be no further consideration of anisotropy resulting 
from the surface energy or wetting potential in this discussion. Further calculations will proceed assuming that the 
surface energy density, 7, nor the wetting potential, W(TL) , depend on V7i or similarly that F s w (H, V7i ) has a purely 
isotropic dependence on V7i. This assumption can be made without affecting any of the qualitative results. 

2.2.1.3 total diffusion potential Having dispensed with the discussion of the various sources of anisotropy, the 
total diffusion potential is stated for the case of 4-fold symmetric elastic anisotropy and a completely isotropic surface 
energy and wetting potential, /ik = /(k, TL) with 

'n[-£ ° (l-£ A sm 2 (26 k ))k + jk 2 + W"(H)] ; case a (Eq. |3)) 
/'( k.H { <-> -£ ° (1 - e A sin 2 (26» k )) k + F^ 2) k 2 + f£ 0) ; case b (Eq. ^) . (23) 
ail — tA sin 2 (2#k)) + bk 2 + c ; case c (general) 
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Table 3: Characteristic wave-numbers, characteristic times and associated dimensionless variables for the three cases 
addressed in Eq. |9]) 





k c 


tc 


a. 





case a 


£o°/i 




k/fc c 


7 W"(H)/£& 


case b 


c /p(02) 
£-0° / " sw 


(f^Y /(DOS*,) 


k/fc c 


p (02) p(20) .pi 
r sw r sw / ^"0° 


case c 


a/b 


b 6 /(vn a A ) 


k/fc c 


do jo 2 



2.2.2 Dynamics 

The dynamics is governed by surface diffusion, just as for the fully isotropic case. It is assumed that the diffusivity 
is isotropic as was done for the surface energy and the wetting energies; thus, all anisotropy in the film evolution 
dynamics comes from elastic effects alone. The possibility and effects of an anisotropic diffusion potential is discussed 
in Appendix [P] (also see [54]). The time dependence of the surface perturbations simply follows Eqs. ( fT3| and ( fB} , 
but with Eq. P3[) used for /(k,7i). As for the isotropic case, appropriate characteristic wave numbers (k c ) and 
time scales (t c ) can be found for each of the three cases along with the associated dimensionless wave vector a and 
dimensionless wetting parameter (3. These are listed in Table [3] The dispersion relation, er^ can be expressed in terms 
of these dimensionless variables (a and 0), giving 

<Jk = o kca = t- l a 2 [a (1 - e A sin 2 (20 k )) - a 2 - 0\ . (24) 

The stability behavior is essentially the same as for the isotropic case with a transition occurring at j3 = 1/4 corre- 
sponding to H = H c - 

2.2.3 Expansion about peaks 

<7k has 4 peaks at (k, 9^) = (ko, ir[n — l]/2) with ko = k c ctQ (Eq. ( fT7| )) and n — 1 ... A. In vector form, there are four 
peaks at 

k„ = k (cos(?r(n - l)/2)i + sin(7r(n - 1) /2)j) . 

Similar to the isotropic case, tJk can be expanded about individual peaks so that in the vicinity of peak n, w a n 
with 

where <r is given by Eq. ( fT8] i, a\\ = a 2 given by Eq. ( fT9l >, and 

<7j^ = 8eACnot c 1 k c . 

In terms of the vector components parallel and perpendicular to k„, kn and k± respectively, 

cr n = CT - 2°"l|( fc ll _ k o) 2 - ^±k 2 ± , 

k\\ = cos[7r(?i — l)/2]k x + sin[7r(n — l)/2]k y , and k± = — sin[7r(n — l)/2]k x + cos[7r(n — l)/2]k y . The time 
evolution of h^in the vicinity of one of the k„ is 

h k (t) sa / lk (0)e t ( CTo -^2(fe||-feo) 2 -i<T_Lfei)_ 

3 Correlation Functions 

Correlation functions and associated constants such as correlation lengths can be very useful for characterizing order. 
In particular, the autocorrelation function (Eq. ( |2~5j )) and its Fourier transform (Eq. ( p6] i) also known as the spectrum 
function can give a very good characterization of dot order (Figs. [6^ and c and[5]3, e and h). The autocorrelation 
function is denoted C A { Ax) where Ax is the difference vector between two points in the x— plane. The spectrum 
function is a function of k, and it is denoted C£. The goal here is to be able to predict these two functions and to 
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describe them quantitatively in a manner that can be used to characterize SAQD order with just a few numbers. The au- 
tocorrelation function is the result of a spatial average over one experiment or one simulation (numerical experiment). 
It is regular and repeatable because it is closely tied to the correlation function and spectrum function that results from 
an ensemble average (Eqs. X and X). These are denoted as C( Ax) and the spectrum Ck respectively. Note that the 
ensemble averaged functions do not have a superscript "A." These ensemble average correlation functions are useful 
in the analysis of stochastic ordinary and partial differential equations. 1331 l56l . From a strictly technical viewpoint, 
the spatial average and the ensemble average are not exactly the same; however, they are closely enough connected 
that it is reasonable to use one as a substitute for the other (Sec. |3.1| and Appendix[3]l. 

In the following, the analysis of SAQD order via autocorrelation and correlation functions is discussed (Sec. [XT}. 



Then, the stochastic initial conditions are discussed (Sec. 3.2 1. Then, the prediction of the Fourier transforms of the 



correlation functions is discussed (Sec. 3.3 I. The real-space correlation functions are presented (Sec. 3.4 1. Finally, 
there are some notes regarding generalizing the analysis method to any dispersion relation that has peaks (Sec. 3.5 1, 
for example, peaks related to broken four-fold symmetry or growth on a miscut substrate. 

3.1 Correlation Functions and SAQD order 

Auto-correlation functions are well-suited for investigating SAQD order. The autocorrelation function is defined as 

^(Ax) = j J dV/i(Ax + x>(x')*. (25) 
Its Fourier transform sometimes called the spectrum ll56ll . spectrum function or power spectrum is 



c * = Wrl rf2Ax ^ !kAxc ( Ax ) = |/lk|2 ' (26) 



where A is the projected area of the film in the x — y— plane. A periodic array of SAQDs leads to a periodic auto- 
correlation function. A nearly periodic array leads to a range-limited periodic auto-correlation function. The ensemble- 
mean of these autocorrelation functions can be calculated, and it is a good predictor of a SAQD order. 

3.1.1 Periodic array 

Consider a perfectly periodic height fluctuation corresponding to a perfect lattice of SAQDs, 

N 

N 



h(x) = ^ ^2 ex P I ik » ' ( x _ x o)] (27) 



plus higher order harmonic, where the dots have a height proportional to ho, N is the degree of symmetry, probably, 
4-fold or 6-fold, xo is a random origin offset. 

, / /2*r(n-l)V . /27r(n-l)\.\ 
k„ = k Q cos —— — '- l + sin — - j , 



N J \ N 

In a linear analysis, the higher order harmonics do not come into play, so they are neglected here. In reciprocal space, 

h N 

n=l 

plus higher order harmonic. The autocorrelation function is found by plugging Eq. ( f27] i into Eq. ( |25| ) and simplifying, 

, x 2 N 

h 



C A (Ax) =1-^1 E exp [ik„ • Ax] (28) 

^ ' 71=1 

plus higher order harmonic. In finding Eq. (|28|), the relation 



J d 2 xV< k '"- k ">' x ' = A5 km k„ = (27r) d S d (k m - k„) 



(29) 
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has been used. <5kk' is the Kronecker Delta, and 5 d (k — k') is the Dirac Delta. Eq. |29) will be helpful whenever it is 
necessary to take an areal average or sum over Dirac Delta functions. In reciprocal space, 

°t - {2] f§2 E ^(k-k m)( 5 2 (k-k„) 

m,n— 1 

h 2 N 

i=l 

plus higher order harmonics, where <5 d (k — k„) = (A/(27r) d )i5kk„ Thus, the order of the SAQD lattice manifests 
itself as periodic functions in real-space (Eq. (|28|) and sharp peaks in reciprocal space (Eq. ([30])). 

3.1.2 Nearly Periodic array 

A nearly periodic arrays shows deviation from perfect order. This deviation is shows itself by broadening of the peaks 
of the spectrum function, C^, and by range limited periodicity of the real-space autocorrelation function, C A (Ax). 
These two measure of disorder are naturally related. 

The disorder in lateral dot size A s ; ze and spacing, A spac i ng are related to each other and to the broadening of the 
peaks in C£ (Fig. [6]a and c). Prior to ripening, the size and spatial order should be related, as the volume of a dot 
should be proportional to the amount of nearby material. If the SAQDs have nearly uniform size and spacing (peak- 
to-peak distance) Lq, the reciprocal space autocorrelation function will be tightly clustered around the wavenumber 
characterizing the dot spacing fco = 2tt/Lq. There are a number of such peaks depending on the system symmetry 
(Fig.[6]a and c), but consider just one. Since the order is not perfect, the peak will have a finite width. Consequently, 
there will be a scatter in the dot size. Since Lq = 2tt /fco, the scatter in dot spacing (A spacnl g) is related to the scatter in 
Fourier components (A^). Taking the derivative of the spacing-wavenumber relation and rearranging, 

^spacing ^ Afc 

L fco 

It is reasonable to expect that the fractional disorder in size (A s ; ze /L s i ze ) is given by a similar (if not exactly the same) 
number. 

Another way to view spatial order (periodicity) is not by dot-dot distances, but the distance over which the dot array 
can be considered periodic. This limited periodicity is evident in the film height autocorrelation function (Eq. |25) and 
Figs. [5]b, e and h). Consider two distant dots. Their position will be completely uncorrelated, so it will be completely 
random as to whether one position corresponds to a peak or a valley. Thus, for a large differences in position the 
autocorrelation function tends to zero. 

C A (Ax large ) = 

Similarly, the mean-square fluctuation of the film height can be large so that 

C A (Ax = 0) > 0. 

The distance over which the autocorrelation function, C" 4 (Ax) decays to is the correlation length, L cor . Thus, L cor 
is a reasonable measure of spatial order. 

The two measures of order A spac j ng and L cor are intrinsically linked. The well known rule of Fourier transforms 
states that the product of the real-space and reciprocal space widths must be greater than or equal to unity. Thus, 
AfeL cor > 1, or A spac i ng > 27tLq/ L cor . Similarly, one can expect that A s i ze ~ L| ze /i cor . Thus, assuming that dot size 
is governed by the amount of nearby material, small dispersions in dot size are only possible if there is long correlation 
length. 

3.1.3 Ensemble Correlation Functions / ergodicity 

SAQDs are seeded by random fluctuations. Consequently, each experiment or simulation must be treated as just 
one possible realization, and the autocorrelation function will be different for each realization. Thus, for analytic 

4 Eq. (29) has been used to help with summation 
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predictions, one must rely on ensemble averages. In G4l . it was assumed that the ensemble average correlation 
function was a good description of a SAQD order, an assumption that was born out by numerical calculations. Now, 
this relation is put on a more solid ground. In particular, it is found that the ensemble correlation functions provide 
good estimates of the auto correlation function and spectrum function produced by any particular realization. First, 
it is shown that the mean value of the film-height fluctuation is zero. Then the method to calculate the ensemble- 
averaged autocorrelation function and spectrum function is presented. Additional mathematical details are presented 
in Appendix [3] 

3.1.3.1 Mean fluctuation It is fairly straightforward to show that the ensemble mean film-height fluctuation is 
zero. The governing dynamics (Eq. ( fT2| i) is invariant upon the substitution h(x, t) — > — /i(x, t). Thus, assuming that 
one does not bias the initial conditions the mean fluctuations must be zero for all time, 

(h(x,t)) = (-/i(x,t)) = 0,and (h k (t)) = 0. 

This is a common situation, and it is most appropriate to characterize the film height fluctuations using the two-point 
correlation function (or simply "the correlation function"). Il55l 



3.1.3.2 Correlation Function The autocorrelation function can be estimated by its ensemble average. Further- 
more, this ensemble average is equivalent to the correlation function that can be easily calculated analytically. These 
relations are first discussed for the real-space correlation functions and then their Fourier transforms. First, the statis- 
tical properties of the autocorrelation function are discussed. Then the statistical properties of the spectrum function. 
Finally, the method to The main results are reported here, and details of derivations are reported in Appendix [E] 

First it is noted that the autocorrelation function averaged over all realizations is equal to the ensemble correlation 
function. 

(C A (Ax)) = C(Ax), where C(Ax) = (ft(Ax)ft(0)) , (31) 

where (...) indicate an ensemble average. Eq. ( [3T] > assumes that the model of film-growth is translationally invariant^] 
This relationship is fortunate, in that it allows one to predict the "typical" autocorrelation function using analytic tools 
that apply only to ensemble averages. 

Second, it is noted that as the area that is used to calculate the autocorrelation function becomes large, the autocor- 
relation function tends towards it mean value, 

C A (Ax)«C(Ax) + 0[A- 1 / 2 ], (32) 

where 0[A -1 / 2 ] indicates statistical fluctuations about the mean value that become smaller and smaller as the area in 
an experiment or the simulation area in a numerical experiment becomes larger. These fluctuations or noise die out as 
A 1 / 2 . For example, the autocorrelation functions in Figs.|5]e and h are very close to the ensemble average autocorre- 
lation functions Figs.[5]f and h, but have random fluctuations that are most visible far from the origin. This property, 
that averaging over a parameter such as position is equivalent to averaging over all realizations, is known as ergodicity. 
Individual realizations are tightly distributed about a "typical" behavior. This tight distribution lends credibility to the 
notion that one can have representative experiments or simulations. Unfortunately, the "demonstration" of Eq. ( |32] > in 
Appendixes not as general as one might like. Rigorously, it applies when the Fourier components of film height (/i k ) 
are independent and normally distributed; however, it is reasonable to conjecture that a relationship like Eq. ([32) holds 
whenever the statistical distribution of film heights is suitably bounded as the boundedness of plays an important 
role in the derivations. 

In reciprocal space, one finds that the ensemble-mean spectrum function is 

(<# ) = C k) (33) 
where C k is defined as the prefactor appearing in the reciprocal-space two-point correlation function. 

C kk = (h k ht) = C\6 d (k - k') = C7 k — ^<W, (34) 



5 A quick survey of literature will find that, virtually all published continuum models of SAQD formation are translationally invariant. 
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where Eq. ( p9| ) has been used. This form for the two-point correlation function in reciprocal space occurs if and only 
if the system is translationally invariant. Eq. ( |3~3"| > is valuable because one can solve for Ck analytically in the linear 
case or using various analytic approximations in the non-linear case. Unlike the autocorrelation function, the spectrum 
function fluctuates greatly about its mean. In fact, the fluctuations are about 100% (Appendi x \E.2\ . These large 
fluctuations result in the commonly observed speckle pattern for the spectrum function C^5 (Figs.|6|a and c). Contrast 
this pattern with ensemble-mean spectrum function Ck shown in Figs.|6]b and d. These speckles can be removed by a 
smoothing operation, and a relation similar to Eq. ([32} results (Appendix |E. 2. 2) . Finally, it s hould be noted that just 
as C^ is the Fourier transform of C (Ax), Ck is the Fourier transform of C(Ax) (Appendix 



E.l 



3.2 Stochastic Initial Conditions 

To model or simulate the formation of SAQDs, it is absolutely essential to include some sort of stochastic effect. An 
initially fiat film /i(x, 0) = is in unstable equilibrium. Thus, to seed the formation of quantum dots, it is necessary 
to perturb the fiat surface. The simplest method to do this is to use stochastic initial conditions with deterministic 
evolution. One can tenuously suppose that white noise initial conditions do not "bias" the ultimate evolution of the 
film. l[57ll Thus, the initial conditions are taken from an ensemble with zero mean, 



(/j(x,0)) = 0. (35) 

and a spatial correlation function, 

C(x, x', 0) = (h(x, Q)/i(x', 0)*) = A 2 6 d (x - x') , (36) 

where the brackets (...) indicate an ensemble average, A is the noise amplitude, and S d (x) is the d— dimensional 
Dirac Delta function. White noise conditions have an infinite amplitude which is not physical. Thus, a minimum 
modification can be made to "cut off" the infinite fluctuations. 



A 2 /_(x-x') 2 
(2 7 r&^ eXP r ^ 



C(x,x',0)= /0 _, 2w/2 exp(- 1 ' (37) 



In the limit bo — > 0, this correlation function reverts to the white noise correlation functions. 
In reciprocal space, 



Ckk'(o) = (Mo)^'(o)) 

= (2n)- 2d J d d x J d d x'e(- lk x+ik ' x ')c(x,x',0) 



^ e-^V(k-k') 



(27r) d 

Letting bo — > 0, the white noise reciprocal space correlation function is obtained. Thus, the initial spectrum function 
is 

The atomic-scale has a small and short-lived influence on the final film morphology (Appendix |F}, but the cutoff 
procedure is useful for choosing a reasonable value of A 2 . It seems reasonable to choose A 2 so that the initial r.m.s. 
fluctuation y/C{0, 0) = (h(0, 0)h(0, O)*) 1 ^ is one monolayer (1 ML). Also, choosing bo = 1 ML as the atomic scale 
cutoff is 

A 2 = (27r) d / 2 (l ML) 2+d , (38) 

where the natural unit 1 ML is, of course, material dependent. 

Using stochastic initial conditions, one can integrate individual initial conditions to obtain representative samples 
and then average over many realizations, the Monte Carlo approach, or one can calculate analytically, the statistical 
measures of the ensemble. The ensemble statistical measures are strongly related to the statistical measures of order 
for an individual realization, so the second approach is opted for here. Thus, the predicted SAQD order is ultimately 
stated in terms of ensemble correlation functions. 
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3.3 Reciprocal Space Correlation Functions 

The reciprocal space correlation function, Ckks and spectrum function, Ck, are calculated for the ID and 2D isotropic 
case and then for the 2D anisotropic case. Generally Ck includes the length scales introduced in Sec. |2. 1.3] as well as 
the atomic scale cutoff b . 

Ock< = (h k (t)hi(t)) = e ( CTk +^'> t <Mo)Mor) 

A * re (^+- k ')M^2(k-k'). (39) 



(27T): 

Without much error, bo can be neglected in the exponential (Appendix IF]). Using Eq. p4| i, the spectrum function is 
then identified as 

C k = ^e***. (40) 
Ck is now calculated for each model: ID isotropic, 2D isotropic and 2D anisotropic. 

3.3.1 one-dimensional 

The one dimensional surface is the simplest, so it is treated first. The spectrum function is simply 

A 2 

r _ ^ r 2a t-^(2a 2 t)(k~k ) 2 

k ~ 2ir 

Ck has a peak at k = ±koi. One can easily read off the correlation length as 

icor = V^2~t = A:-V2(3a -4/3)(t/t c ). (41) 

so that 

A 2 

n. - r 2a t-^Ll r (k-ko) 2 

This approximation is valid when koL cor ^> 1. In terms of k x , 

3.3.2 2D isotropic 

The 2D isotropic case is very similar; 

A 2 

r< — p 2cr ot-5iL( fe - fe o) 2 (aj\ 

Ok- (2?r)2 e , (42) 

where L cor is the same as in Eq. ( |4T| i. It has maximum that forms a ring in the k— plane as graphed in Fig.[6]b. 

3.3.3 anisotropic 

The anisotropic spectrum function is 

C k = ^^e^o^e-i^ii-M 2 -^*!, (43) 

v y n— 1 

where 

L|, = = fc-V(6ao-8/3)(Vi c ), (44) 

L_l = \/zVl< = fc-Vl6eao(i/i c ), (45) 

= cos[7r(n — l)/2]fc x + sin[7r(n — l)/2]fcj,, and A;^ = — sin[7r(n — l)/2]k x + cos[7r(n — 1) /2]k y and it is graphed 
in Fig.|6]d. This approximation is valid when k L\\ ^> 1 and k L± ^> 1. 



19 




-0.3 -0.2 -0J 0.1 0.2 0.3 -0.3 -0.2 -0.1 0.1 0.2 0.3 

kj -dim" 1 ) k x (ma~ l ) 



Figure 6: C£ and C^for Ge/Si as discussed in Sec. [4] (a,b) 2D isotropic surface. Eq. ( |42] i is used for Ck. (c,d) 2D 
anisotropic surface. Eq. ( |43] l is used for Ck- 
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Figure 7: ID isotropic surface in real space for Ge/Si as discussed in Sec. [4] (a) Example of h(x) plotted over a length 
of 8L cor . (b) corresponding reals space correlation functions plotted for range ±4L cor . Filled plot is an example of 
C A (Ax). Solid UneisC(Ax) (Eq. |46}). 



Loosely speaking, one can argue that the isotropic case is similar to letting ca — * in Eq. ( |45] l so that the 
perpendicular correlation length is always regardless of time. A more conservative approach would be to argue that 
Lj_ w 2n/k for the isotropic model via inspection of Figs. |6ja) and (b). Even still, the more conservative result 
guarantees that the perpendicular correlation length will always be the same as the dot spacing; thus, it will always 
limit SAQD order to the first nearest neighbor at best. 

3.4 Real Space Correlation Functions 

The real space correlation functions C(Ax) are now calculated for the ID and 2D isotropic cases and the 2D elastically 
anisotropic case. 

3.4.1 one-dimensional 

In one dimension, 



/oo 
dk x e ikxAx Q 
-OO 



DC 

-V L e 2<Tot -5A* 2 /£l2cos(fc x). (46) 



Thus, C(Ax) has a damped periodicity indicating that it is imperfectly periodic (Fig. [7]) 

3.4.2 2D isotropic 

In two dimensions with elastic isotropy, 



C(Ax) = J d 2 ke ,;k Ax C k 



(2^ ' ^ k 



o Jo 



Performing the angular integration first, 

,2 



iL 2 (k-k ) 2 



C(Ax) = ^-e 2<Tot J dk kJ (kAx)e-^ 
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where Jo is the zeroth Bessel function. In general, this integral is best performed numerically; however, it can be 
solved in two important cases: Ax — * and L cor — * oo (corresponding to long times). In the first case, 



C(Ax) = — e 2CTot / dkke-^ k - k ^. 
2tt Jo 



Under the same conditions that Eq. ( pi) is valid (koL cor 3> 1), the lower limit of the integral can be approximated as 
—oo so that 

C(Ax = 0) = ^! fc ° e 2CTot . (47) 



2ttL, 



cor 



\L 2 „(k-k ) 2 



This function gives the mean square surface height fluctuation. In the second case where L cor — ► oo, e 2 
(27T) 1 / 2 L~lS(k - ko), so that 

C(Ax) = e 2CTpt J (fcoAx) . (48) 

V 2-KLccr 



This correlation function is the most ordered case for a 2D isotropic surface. It is graphed in Fig. [5};. 
3.4.3 anisotropic 

To find the real-space correlation function for the elastically anisotropic case, it is best to find the contribution from 
each peak and then sum so that 

A 2 4 

C(Ax) = 7^2 e2CT0 * E C"( Ax ) (49) 

^ n ' n=l 

where 



C"(Ax) = y d 2 ke' k ' x e _ 5 i ii(' C|l_feo ) 2_ 3 i i fc i. 



Ax can be decomposed into the directions parallel and perpendicular to k„, so that Ax\\ = cos(7r(n — l)/2)Ax 
sin(7r(n - 1) /2) Ay and An = - sin(7r(n - l)/2) Ax + cos(7r(n - 1) /2)Ay. Thus, 



C"(Ax) 



dk\\ e lk « Ax «-^ L «( k ^ ko ) 2 ^j dk ± e 



ik\Ax\—kL 2 \ k 2 . 



277 - e -H x V L l+ x ±/ L2 L) e ikoX " . 



L\\L\ 



Plugging into Eq. ( |49) , 

A 2 



C(Ax) = ^_ e 2-o* r e -i(x a /iS+» 2 /Ll) cos (fc oS ) + e-K^/ii+j/Vif) cos (fc 2/) 
ttLi\L± 



(50) 



3.5 Generalizability 

The dynamics and analysis used here were for a specific model, but the general procedure for analyzing the order 
resulting from a linearized model should hold for any model with well-separated peaks in the dispersion relation, a^. 
The procedure to follow is: 

1 . Generate the dispersion relation, as some function of k. 

2. Find the peaks in the dispersion relation, k„, (n = 1 . . . N) 

3. Expand about the peaks to generate the peak values, a n , and local Hessian matrix, 

d 2 



ij dkidkj 



Ok 



k=k„ 
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The spectrum function is then approximately 

A 2 N 

C k (i) 



,2o-„i 



exp 



t (k-k„)-# n .(k-k„) 



(51) 



4. Find the Eigenvalues of the local Hessian matrix, (H n )j and (H n ) u . They should be negative, if there is a 
peak at k„ 

5. Use the eigenvalues to determine the correlation lengths, (L n )j = v/2 |(i? n ) J | i and (L n )jj = \/2 £. 
The real-space correlation function is 



A' 



(2tt) ^ 4V(^«)/ (S3 



.e 2CT "*exp 



x • H n 1 ■ x 
4t 



(52) 



The "goodness" of these approximate forms requires that (L n )j and (L n )Jj be much less than the spacing 
between peaks in the correlation function so that the gaussians do not overlap greatly. A reasonable test for 
this no-overlap condition is ||k n || (L n )j -C 1 and ||k n || (in)// "C 1, assuming that the peaks are not large in 
number or very closely spaced. 



4 Order Predictions 

The real-space correlation function formulas (Eqs. ( |46l ), |47|, and |50]l) and correlation length formulas (Eqs. ( [41] ), ( |4~4~t > 
and ( f45] >) can now be used to estimate the order of SAQDs. Ge on Si is chosen for this example because this system has 
received the most attention from theoretical work [58 38 31, 18, 39||4T][22]|25]|26] and others], and it is the simplest 
since it involves the diffusion of a single species. The procedure described below tries to predict the amount of order 
when an initial atomic-scale fluctuation becomes "large". "Large" is taken to be greater than atomic-scale. Beyond 
this point, one would expect non-linear terms to become important. An example is presented for Ge on Si at 600K to 
compare and contrast the 2D anisotropic results with the ID isotropic and 2D isotropic results. The predictions are 
also compared with a linear numerical calculation on a discrete reciprocal -space grid to test the approximations made 
and to illustrate the relation between the surface profile (/i(x)), the example autocorrelation functions (C A (x) and 
Cj^ ) and the ensemble correlation functions (C(Ax) and Ck). Figs.|6]|7]and|5]show these results. Finally, the relation 
between average film height and order is investigated. 



4.1 Ge at 600K 

The formulations for the three discussed cases are implemented for Ge/Si at 600K. The correlation lengths are esti- 
mated for the end of the linear regime where fluctuations become large (greater than atomic scale). First, appropriate 
physical constants are used to give the corresponding correlation length and correlation functions vs. time. These in- 
clude an initial average film height Ti and a white noise amplitude A (Eq. ([38|). These initial conditions approximate 
a film at the beginning of an anneal that immediately follows a rapid deposition. The time t\ mge is found by solving for 
the time where the mean-square fluctuations are atomic scale, (/i(x, t) 2 ) = C(Ax = 0) = 1 ML 2 . At this point, the 
correlation lengths are calculated. 

Physical constants for the 2D anisotropic calculation are taken as follows. The elastic constants for Ge at 600 K 
are c n = 1-199 x 10 12 , c 12 = 4.01 x 10 n (from C : s = 3.991), c 44 = 6.73. ED Using a Ge = 0.5658nm and 
a S ; = 0.5431nm, it is found that e m = 0.0418. Using the procedure from (Appendix^, M = 1.332 x 10 12 dyn/cm 2 . 
£o° = 4.96 x 10 9 erg/cm 3 , and £450 = 4.35 x 10 9 erg/cm 3 , giving = 0.1236. The atomic volume is fl = 
2.27 x 10~ 23 cm 3 . The estimated surface energy density is 7 = 1927 erg/cm 2 . The wetting potential is estimated 
by picking a plausible critical surface height, TL C « 4 ML = 1.132 nm and setting W(H) = SqoHI / (8jH) = 
2.315 x 10~ 6 /7Y erg/cm 2 . The resulting characteristic wave number is k c — 0.257 nm -1 . The initial film height is 
taken to be H = H c + 0.25 ML = 1.203 nm and then allowed to evolve naturally. Thus, (3 = 0.208, a = 0.5658, 
k = 0.1456 nm" 1 , a = 0.1192/i c , a\\ = 0.864/(fc 2 i c ), a ± = 0.559/(fc 2 i c ), L\\ = 0.744fc ( 7 1 (t/t c ) 1 / 2 , and 
L± = 0.599/c ~ 1 (/;/i c ) 1 / 2 . The unspecified diffusivity has been absorbed into the characteristic time t c . From Eq. ( f5~8] >, 
A 2 = 0.0403 nm 4 , and Eq. Q gives 

C(0) = (1.223 x 10- 3 i c /i) e ' 023854 ^ nm 2 . 
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The initial infinitely rough surface undergoes a smoothing described by the t c /t factor. Then the surface roughens due 
to the exponential. The initial divergent roughness is an artifact of the non-physical white noise with the atomic scale 
cutoff &o neglected (Appendix [F]). The time for the fluctuations to become "large" again are found by setting 

C(0) = hl ge (53) 

where /iiarge = 1 ML = 0.283 nm. The solutions are t\ = 0.01527t c or t% — 430i c . The first solution is discarded 
since it is due to the non-physical white noise. At t\ Mge = t 2 , Ln = 105.8 nm, and L± = 85.2 nm. Taking L± as 
more limiting, the correlation spans about n = k^L^/n = 3.95 islands across. The corresponding reciprocal space 
(Eq. ( |4"3"j )) and real-space correlation function (Eq. ( |50] >) are shown in Figs.|6]d and[5]f respectively. 

A corresponding numerical experiment is performed. A periodic surface of size I = 96(27r/fco) is used. Random 
initial conditions consistent with Eq. ( [38] l are used for k— space points on a square grid bounded by k x , k y — ±2fco- 
The relation between discrete and continuous Fourier components is used, (/ik) discrete = [(27r) d /A]/ik- Eqs. ( fT3j ) 
and ( fl4] i are used without any additional approximation to find at time t — tiarge- The resulting C£, a portion of 
the height profile /i(x) and C A (Ax) are plotted in Figs.[6jc),[5jd) and[5jf) respectively. 

Similar calculations can be performed for the one-dimensional and two-dimensional elastically isotropic cases. 
Isotropic values used previously [58, 24| are about E = 1.361 x 10 12 dyn/cm 2 and v = 0.198 giving M = E/(l — 
v) = 1.697 x 10 12 dyn/cm 2 and E = 2M(1 + v) = 7.10 x 10 9 erg/cm 3 . Using the same critical surface height, 
He = 4 ML, W(H) = 4.74 x 10~ 6 /?i erg/cm 2 . The resulting characteristic wave number is k c = 0.368 nm -1 . If the 
film is grown to H = H c + 0.25 ML = 1.203 nm and then allowed to evolve naturally, [3 = 0.208; thus, a a = 0.5658, 
k = 0.208 nm -1 , a = 0.1192/t c , a 2 = 0.864/(fc 2 i c ), and L ml = 0.744A: r 7 1 (i/t c ) 1 / 2 . In one dimension, Eq. (|46]> 
is used to find the mean square height fluctuation. Using Eq. ( f3~8| ) with d = 1, A 2 = 0.0568 nm 3 , and 

C(0,i) = 0.01271(Vi c r 1/2 e a0238t/tc . 

Setting C(0,t) = (1 ML) 2 = 0.0801 nm 2 , t x = 0.0252i c , and t 2 = 186.9t c . At t 2 , L cor = 48.8 nm, and n = 
koL COI /TT = 3.24, so about 3 dots in a row should be well correlated. The corresponding numerical calculation of 
size I = 96(27r/fco) is performed. A portion of h(x), C j4 (Ax) and C(Ax) are shown in Fig. [7] In two dimensions, 
Eq. (|47]> is used to find (h(x, t) 2 ), 

C(0,t) - 9.40 x lO-^t/Q-^e - 0238 ^. 

Setting C(0, t) = 0.0801 nm 2 , t x = 1.376 x 10~ 4 i c , and t 2 = 306t c . At i 2 , L COI = 62 .4 nm, and n = k L cor /Tr = 
4.14, and correlation is expected to extend about 4 dots. However, it should be noted that this correlation is not 
lattice-like. Corresponding numerical results and ensemble correlation functions are shown in Figs.[6]and[5]a-c. 

4.2 General case of (3 

In 1 24 1 it was suggested that allowing the film to evolve with /3 close to the stability threshold could enhance the SAQD 
correlation. It is interesting to note what happens for different values of /3. Similar analytic and numerical calculations 
are performed for the large film-height limit, (3 = 0, for the 2D anisotropic Ge/Si surface. For (3 = 0, <i arg e = 40.3i c , 
L± = 30.0 nm, and n = k^L^jn = 1.84, so one to two dots in a row are expected to be well correlated. /i(x) 
and real-space correlation functions are shown in Figs.[5^-i. The range of order is significantly less than for the case 
(3 = 0.208 (Sec. |4. l} . For Si/Ge at 600K, the 2D anisotropic predictions for i^ge and L± are shown in Fig. [8] In 
general, the closer j3 is to the critical value 0.25, the longer the correlation length. One can manipulate equation p3[) 
to find that t\- drgs: /t c varies approximately but not exactly as ((3 — 1/4) -1 x ln[/i 1 2 argeV / eA/(A 2 fc 2 )]. Consequently, 
Lj_ ~ (/3 — 1/4)~ 1//2 . Furthermore, the appearance of /ii mg e and A 2 inside the logarithm shows that the final order 
estimates are not overly sensitive to the guesses for A 2 and hf . The divergence of L± with (3 — 1/4 is initially 
encouraging, but it is clear that for the parameters used for Ge/Si, subatomic control of the film height is needed to 
yield significantly enhanced long range correlations. Also as one approaches this threshold, one can probably expect 
thermal activation to nucleate subcritical SAQDs whose effect on supercritically formed SAQDs is uncertain. There 
should be some interesting phenomena at the the Ti — * Ti c . 
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t lai e and L ± vs. (3 (Ge/Si) 




0.05 0.1 0.15 0.2 0.25 
P 

Figure 8: ti arge and L± vs. (3 for Si/Ge using the 2D anisotropic model as described in Sec. [4] Units are normalize to 
characteristic time t c and predicted number of correlated dots (n = koL±/ir). 

5 Discussion/Conclusions 

The order of epitaxial self-assembled quantum dots during initial stages of growth has been studied using a common 
model of surface diffusion with stochastic initial conditions. It has been shown that correlation functions of small 
surface height fluctuations can be predicted analytically using corresponding ensemble average correlation functions. 
These correlation functions are characterized by correlation lengths that can be predicted by analytic formulas given 
certain reasonable assumptions about the diffusion potential and the height and lateral scale of initial atomic scale 
random fluctuations. Thus, the linear model of film surface height evolution via surface diffusion has enabled analytic 
predictions of epitaxial SAQD order that are valid for small film height fluctuations. To what extent the initial degree 
of order persists into later stages of growth remains to be studied, but the order of initial stages should certainly have 
a strong influence on final outcomes. Furthermore, the linear analysis should provide insight into the less tractable 
non-linear behavior. These predictions of SAQD order have been used to investigate the role of crystal anisotropy and 
initial film height. 

Crystal anisotropy has been shown to play an important role in enhancing SAQD order as observed in previous 
numerical simulations continuum and atomistic numerical simulations. B31 l37l l44l l45l If a four-fold symmetry is 
assumed for the governing dynamics, the effect of crystal anisotropy to linear order is felt through elastic anisotropy 
alone. It is shown that elastic anisotropy is required to produce a lattice-like structure of SAQDs. The enhanced spatial 
order should in turn lead to enhanced size order, a consequence that must be confirmed with non-linear studies, but 
appears to be true based on the present available literature. 

The role of initial film height has been shown to greatly influence order. Growth near the critical film height for 
dot formation can enhance order. This order enhancement comes from increasing the duration of the linear small- 
fluctuation stage of growth. In fact, the predicted correlation lengths diverge when the initial film height approaches 
the critical film height from above. Achieving large correlation lengths in this manor is of course practically limited 
by ability to control film heights to subatomic accuracy. Additionally, one should be careful when interpreting the 
continuum model in such a context, as the effect of atomic discreteness might be greater at the transition film height. 
Finally, it is likely that additional randomizing effects of thermal activation will effectively cut off this divergence 
when the critical film height is approached from below during deposition. 

Finally, the presented method may be useful as a first step in the analysis of methods to enhance SAQD order. It is 
reasonable to suppose that under some circumstances initial growth stages will be very important while for others they 
will not. For example, prior work on vertical stacking appears to confirm the presented ordering mechanism. |44|. 
Vertical stacking not only achieves vertical correlation of dots, but each layer is more ordered horizontally than the 
one below. Additionally, a "growth window" was found, whereby to achieve enhanced order, the evolution of each 
layer be terminated before ripening begins. The reported simulation [44] supports the following scenario for SAQD 
order development. Order is enhanced during the small fluctuation stage as described here. Once the fluctuations are 
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sufficiently large, the seeded dots evolve towards their equilibrium shapes. Finally, the dots begin to ripen and order 
diminishes. Order is transfered via strain to the next layer so that the next layer gets a head start on its initial ordering. 
Thus, the multiple layers of dots effectively draws out the linear growth stage. It may be possible to modify the present 
model to predict the correlation length of each SAQD layer. 



A Diffusion Potential 

The diffusion potential is calculated in terms of the film height H that is a function of the in plane coordinates x = 
xi + yj. The elastic and surface energy portions of the diffusion potential can be found in ifTSl 

Melast(x) = £M X )> and Vsurf = -^7 K ( X ), 

where il is the atomic volume, w(x) is the elastic energy density at the film surface, 7 is the surface energy density, 
and k is the total surface curvature. However, other calculations need to be included: 

1 . /i wet for the two wetting potential cases, Eq. |3) and d3), 

2. and /.t sur f and /i wet when the surface energy density 7 and wetting energy density W also depend on surface 
orientation. 

Before these case are addressed, a general form for the diffusion potential is justified. 



A.1 General Form \i = fi£F/5ft(x) 

The diffusion potential, /x(x), is the change in free energy, T, when a particle is added at a position, x. Note that /x(x) 
and T are relative energies. They can be used to compare the binding energy of one site on the surface in comparison 
with another site, but should not be interpreted as an absolute binding energy or total formation energy of the surface. 
If a particle has a volume f2, then the diffusion potential at x is related to the variation of free energy with volume, 

5T= Q' 1 J d d xfj,{x)5V(x), (54) 

where 5V(x) is the volume variation at x. Calculating SV(x), V — J d d xH(x). Therefore, 5V(x) = STC(x). 
Substituting into 6f (Eq. ((54)), ST = fT 1 J d d xn(x)5H(x) or /i(x) = nST/SH(x). 



A.2 Simple Model 

Starting from Eq. ((2), /i(x) is found by taking the variational derivative, 

S f 

Meiast.(x) = — - / d d xdzu{H](x, z) = flu; (x) 

Ort(X) „/ vo i um e 

where the "[%]" indicates that the elastic energy, oj, is a nonlocal functional of the film height Ti, and ui(x) = 
uj[H] (x, H(x)), the elastic energy density evaluated above lateral position x at the free surface (z = H(x)). See lfT5ll 
for details of the derivation. The surface energy diffusion potential is 

/W(x,i) = il-^—Jdfix [l + (VH(x)) 2 ] 1/2 7 

= -OV-[1 + (VH(x)) 2 ] 1/2 7 = -0 7 k(x). 

The wetting energy diffusion potential is 

Mwet(x) = fi^y J d d xW(H(x)) 

= nW(H(x)) 
Putting these three terms together, one obtains Eq. ((3) 
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A.3 General Model 

Consider the general form for the combined surface energy and wetting potential, 



F sw = J AF SM (H(x),VH(x)) 



as in Eq. Q so that the free energy is an integral over the x— plane of an energy density that depends on 7i(x) and 
VH(x) locally. The corresponding diffusion potential is 

M(x) = tt^g = n [^£ 0) (H(x), VH(x)) - V • F&°>(W(x), VH(x)) 

B Linearized Diffusion Potential and Anisotropy 

The linearized diffusion potential ^ij n> k is found by finding n(x) to first order in height fluctuations (h), to get /iii„(x) 
and then taking the Fourier transform to get /iu ni k- The linearization of the simple isotropic diffusion potential corre- 
sponding to Eqs. |2) and Q was discussed in Sec. 2.1.1.1 Here, the more general diffusion potential corresponding 



to Eqs Q and |5]) is linearized and then applied to the anisotropic simple model and the anisotropic general model 
Only the surface and wetting parts of the diffusion potential are discussed in this appendix. See ref. fl5l . Sec. 2.2.1.1 
and Appendix |C| for discussion of /i e iast.- 

B.l Linearizing the simple model 

Consider a wetting potential and diffusion potential that both depend on the film height gradient V7i, 7 — > 7(V7i) 
and W(TL) — ► W(H, V7i). Starting from Eq. (|6]l and expanding to second order in the film height fluctuation using 
W(x) = H + /l(x) (Eq. 0), 

1 + (VW) 2 ] V2 \(VH) = (l-^(Vhf + ...^j (7 + 7' • V/i + 7 " : VhVh+...) 

= 7 + 7'- Vh- ^7(V/i) 2 + 7" : VhVh + 0[h 3 } 
where 7 is 7(0), and the primes indicate the derivatives with respect to the surface height gradient. 

i = 9vk7(VH)| vh=0 , and 7" = d VH dvHj(VH)\ vn=0 . 

Taking the derivative with respect to V/i results in a tensor of rank equal to the order of the derivative because V/i is 
a vector (tank 1 tensor). Taking the variational derivative, /x sur f.(x) = £7(5jF surf ./<5/i(x), 

/W.,iin(x) = n [ 7 V 2 />(x) - 7" : VV/i(x)] 

The term with 7' vanishes because it is the divergence of a constant (V • 7')- Taking the inverse Fourier transform, 

AW., iin,k = Q {-jk 2 + k • 7" ■ k) Ai k . (55) 

The first term is isotropic. The second term is parameterized by a rank 2 symmetric tensor. 

Going through the same process, one finds essentially the same result for an orientation dependent wetting energy. 
The step details are so close to the details for linearizing the more general form, F SW (TL : V7i), they are deferred to 
(Appendix |B.2| i. One finds that 

(V (20) + k • W (02) • k) . (56) 

where W {mn ^ = d%d^ H W(H, VH)\ H=n VH=0 is the m th and n th derivative of the wetting energy density with 
respect to TC and V7i evaluated for a perfectly flat film of height 71. W^ mn ^ is a tensor of rank n. 
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B.l.l isotropic case 

In the isotropic case, 7" — * 7"!, where I is the identity operator, and 7" is a scalar. Similarly, W^ 02 ) — > W^ 02 'I. One 
thus gets for the combined surface and wetting parts of the diffusion potential, 

M«»,nn,k = n [(-7 + 7" + ^ (02) ) fc 2 + V^ (20) ] fc k . 

Thus, in the isotropic case, the linear order effect of introducing a surface orientation to either the surface energy or 
the wetting potential is simply to change the apparent surface energy density by 7 — > 7 — 7" — W^ 2 \ 

B.1.2 anisotropic case 

The surface and wetting parts of the diffusion potential (Eqs. |55) and ( |56) )) can admit only a limited anisotropy. 
They both contain rank 2 symmetric tensors, 7" and W^ 02 ) in the x— plane. For a two-dimensional surface, this 
means that they can either have two-fold-symmetric (rotations by 180°) anisotropy or none at all. Thus, for the case 
considered in Sec. |2.2.1.2| four-fold-symmetric anisotropy , the surface and wetting parts of the diffusion potential 



must be completely isotropic. As discussed in Sec. 2.2.1.2 the (100) surface of zinc-blend structures, such as the 



mentioned Ge, Si, InAs and GaAs present a rather complicated situation. For simplicity, it is assumed here that the 
surface and wetting energies are at least four-fold symmetric. Consequently, they are completely isotropic. 

Finally, it should be noted that if F sw depends on higher order derivatives, then the discussion is greatly compli- 
cated and a larger class of anisotropic terms is admissible. For example, when F sw — > F sw (H, V7i , WW, W V7i, . 
is expanded about K(x) = Ti to quadratic order in h, it would contains tensors of rank 6 and maybe even higher. 

B.2 Linearizing the general model 



The elastic part of the linearized diffusion potential was discussed in Sec. 2.2.1.1 and Appendix [C] . Eq. |56) can be 



found by using all of the following steps with the substitution F sw — > W. The surface-wetting part of the diffusion po- 
tential /i(x) is found by expanding F sw to second order in the film-height fluctuation, h, and then taking the variational 
derivative. Expanding F sw about h = and V/i = 0, 

F sw (H+h,Vh) = Fj°°) + F&°) h + F(°J) • Vh + hF<£> • Vh . . . 

■■■ + \F^h 2 + \¥ s ^ -.VhVh + 0[h\ 

Note that in this expansion, all the Fsw terms are constant with respect to h and depend implicitly on the average 
film height, H. The first index indicates the m th derivative with respect to h. The second index indicates the n th 
derivative with respect to V/i. The derivatives are evaluated for a perfectly flat surface of height H. Thus, 

FjS" = d%d$ n F sw (H, VH)\ H ^n=o ■ 

Since V/i is a vector in the x— plane, Fsw^ is a tensor of rank n. Taking the variational derivative of T sw — 
J d d x F SW (H, V7i) and keeping terms to order h 1 , 



ST 
6h(x) 



- V • F( 01 > + F^h - V • ( F(° 2 > • Vh) 



Note that the FiS, ' term vanishes because it is constant, and the F^'' term vanishes upon simplification. Additionally, 
the Fsw can be neglected if one enforces the condition that the film-height fluctuations do not add or subtract material 



from the surface, namely that J d d x 5h(x, t) = 0. Alternatively, one can discard it in anticipation of taking the gradient 

of the diffusion potential, since it is a constant. The t 
constant. Multiplying through by the atomic volume, 



of the diffusion potential, since it is a constant. The term V • F,™' > = for the same reasons, or because FsuP is a 



Mlii 



(57) 
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B.2.1 isotropic case 

In the isotropic case, F^™ must be proportional to the identity so that F £2, = Fi" 2 'l; thus, 
Taking the inverse Fourier transform of this equation, 



/^si/J,lin,k ^ 



F (20) , F (02),2 



This gives case b in Eq. Q. 



B.2.2 anisotropic case 

If the surface is anisotropic, then Fi™ in Eq. ( |57| ) is a rank 2 symmetric tensor in the x— plane. Thus, it can have two 
distinct eigenvalues, and automatically has 2-fold rotational symmetry (rotations by 180°). If any other symmetry is 



assumed such as 4-fold symmetry (rotations by 90°), then F 
transform, 



(02) 



must be fully isotropic. Taking the inverse Fourier 



/^su>,lin,k 







F m + k . f(02) . k 

Old 1 O+JJ 



In Eq. ( |23j >, case b, it is assumed that there is four-fold symmetry, resulting in a surface-wetting part of the diffusion 
potential that is completely isotropic. 



C Elastic Anisotropy 

In principal, the anisotropic elastic energy is found in the same fashion as the isotropic elastic energy. [15| The 
fiat film, initially in a state of biaxial stress, is perturbed by a small periodic surface fluctuation of amplitude ho. An 
appropriate elastic field is added to satisfy the perturbed traction-free boundary condition at the free surface. Finally, 
the elastic energy is evaluated at the free surface to first order in ho. The coefficient ho is the sought after o»k- The 
equations themselves are cumbersome and best solved using a numeric implementation, so an abstract procedure for 
calculating cu^ is outlined here, is found for k = 1 but arbitrary 8^. 
Let the surface have a height variation 

fc(x) = h e lkx . 

To first order in ho, the surface normal is 

n(x) = -ikh e ikx i + k. 

The elastic energy needs to be calculated to first order in ho- To find the elastic energy, it is necessary to find the 
perturbing elastic field to first order in ho- 
The initial unperturbed stress state is 



Cm 





" 





Cm 















where <r m = (en + en — 2cf 2 /cn) e m . Note that this stress state is isotropic in the x~ y-plane and thus independent 
of rotations about the vertical axis. Under this stress state, a flat surface is traction-free. With the height perturbation, 
the traction is 

tj = (n • a m ). = -~ikhoMe m 5 n e ikx . (58) 

Next to find the perturbing elastic fields. These are not isotropic in the x — y— plane, and it is necessary to take 
into account the angle. First, the 3x3x3x3 elastic stiffness tensor Cyjy is constructed for the cube orientation from 
the compact 9x9 matrix cy. The tensor representation aids in rotation. The stiffness tensor is then passively rotated 
in the x — y~ plane by an angel 9^, 

3 

Cijkl(0k) = ^ R{0\t)imR{6]i)inR{0\t)kpR{9\t)lqCmnpq 

m,n,p,q—l 
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where 



cos(6>k) si n (^k) 
— sin(0k) cos(^k) 
1 



This passive rotation of Cijki is equivalent to actively rotating the wave vector k = ki by 9^. 

The appropriate form for the perturbing displacement field is found. Assume a displacement of the form 

u i (x,y,z) = U i e k ^+ KZ \ 

where n can have a complex value. The elastic equilibrium equations are 

3 Q Q 



i,k,l=l 



where 



Y, C dh,K)vA k 2 e k( - tx+KZ ^=0 
\i=i J 



Cji(9 k , k) = 2_j c ijki(dk)(iSu + S l3 K)(iS k i + 6fc3«0- 



(59) 



i,k=l 



Factoring out k 2 e k ^ x+KZ ', the part in parenthesis must be identically zero. 

To obtain a non-trivial solution, the determinant of Cji(9 k , ft) to zero. Six complex values of k are found. The 
values of k with Re[«] < are discarded since the corresponding displacements blow up as z — > — oo. Each of 
the remaining values k = k p with p = 1 ... 3 is substituted back into Cji(9k, k), and Eq. ( |59] > is solved to find the 
corresponding eigenvectors, Uf. The total displacement is thus 



ui(x, y, z) = ie m h A p U[e k{lx+KPz \ 
P =i 

where it is assumed that the perturbing elastic displacement field is proportional to h and a m , and the factor of i is 
put in for convenience. The coefficients A p can be found from the traction-free boundary condition at the free surface. 
The traction formula is 



9 

n lCljk i(9 k )—ui(x,y,z) = ike m h ^ niC ljk i{e k )A p Uf{i5 kl + K p 5 kz )e Htx+KPz) (60) 



i,k,l=l i,k,l,p=l 

The traction is already proportional to ho- Thus, all terms in the sum must be kept to zeroth order in ho so that 

h(x) = 0, and n(x) = k. 

Thus, plugging z = to Eq. ( |60l >, 

3 3 



tj = ike m ho Y Z ( ic 3jii(#k) + K p c 3i3 i(^k)) A p Ufe tkx . 



(61) 



=i i=i 



Since the total traction (Eqs. ( |58] > and ( |6*Tj >) must be zero, the coefficients A p are found from 



where 



K JP = Y, (ic 3j ii(^) + « P c 3 , 3i (0 k )) Uf, 
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and 

R 3 = M5jt 

for j = 1 ... 3. It is worth noting that only for the symmetry directions, 6*t = 0° and 0k = 45° is the strain purely 
plane-strain as it is for the elastically isotropic case. 

The elastic energy at the film surface is found to order 0(h a ). If the stress and strain are expanded to first order in 

ho, a = <7o + cti, and e = €q + ei, then 

U = : c : e = ^a : e Q + cr : ei + 0(h%). 

Thus, 

U = U + Me m ((ei)n + (£1)22) 

du 3 
(ei)n = -7^7 = -emkhp A p Uf. 

P =i 

(£1)22 = du 2 /dy = 0. Thus, 

f/ = C/o -£e k kh a e lkx 

where 

3 

P =i 

where A p and J7f are implicitly functions of 9^. This procedure has been used to find the values of £o° an d £45° for 
Tableland Sec. g] 



D Diffusional Anisotropy 

In general, the surface diffusivity can depend on the film height 7i(x) and the surface orientation V7i(x) so that the 
surface current is 

Jg(x) = D(W(x), Vi?(x)) • Vsfi(x) 

where Vs is the surface gradient, and D is a rank 2 tensor in the two-dimensional space tangent to the film surface at 
x. Linearizing the surface current about a fiat surface, 

J s (x) = D(W) • VMim(x) 

where the diffusivity must be evaluated for h — and V/i = 0, since /J,i; n (x) is already proportional to h(x). T he lin 



earized diffusivity is a symmetric rank 2 tensor in the x— plane. Thus, it is similar to F sw discussed in Appendix B.2.2 
It is automatically either two-fold symmetry (rotations by 180°) or it is completely isotropic. In Eq. < |23) , four-fold 
symmetry of the surface is assumed. Thus, the diffusivity must be completely isotropic; D — * T>, a scalar. Sec- 
tion 2.2.1.2 and Appendix B.2.2 contain discussions of the symmetry properties of the various rank 2 tensors that 
appear in the linear evolution equations. A limited case of diffusional anisotropy has been modeled via kinetic Monte 
Carlo technique. l54l 



E Correlation Functions 

E.l Mean Values 

Equations ( |3"T| and ([33) are central to the presented analysis. Here, they are derived. The two-point correlation func- 
tions for a stochastic system are introduced. Then, the average of the autocorrelation function is taken and expressed 
in terms of the two-point correlation functions. Finally, this average is simplified using the translational invariance of 
the system (governing equations and ensemble of initial conditions). 
The two-point real-space space correlation function is 

C(x,x') = </i(x)ft(x')*), 



31 



and the reciprocal space correlation function is 

= (h k ht,) . 

These are related by the double Fourier transform, 

CW = j^m J d d xd d x'e- 4k x + lk ' x 'C(x,x'); (62) 
C(x,x') = /d d kd d k'e ikx - lk ' x 'C kk /. (63) 

These ensemble correlation functions can be used to give the ensemble-mean autocorrelation function and spec- 
trum function. In real space, 



(C A (Ax)) = i /" dV (/i(Ax + x')/i(x')) 
= j /rf 2 x'C(Ax + x',x'). 



(64) 



(Q 4 ) = ^ («> = ^C kk . (65) 

Fortunately, the translational invariance of the system simplifies these relations. Inspecting the governing equations 
and invoking the translational invariance of the stochastic initial conditions, the resulting ensemble and its statistical 
measures must also be translationally invariant. Thus under the translation by x', 

C(Ax + x',x / ) = C(Ax,0) = C(Ax), (66) 

so that the independent variable is reduced to just the difference vector Ax = x — x'. This relation can be used to 
simplify both the real and reciprocal space relations. 

The real space relation simplifies as follows. Inserting Eq. d66) into Eq. d64), 



(C A (Ax)) = j J d 2 x' C(Ax, 0) = C(Ax). (67) 
The reciprocal space relation (Eq. ( |62| )) simplifies to 

CW = C*6 2 (k - k') = C k -^<W, (68) 



where 



' ,.l y! j rf 2 Axe- lk ^C(Ax). 
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One can see immediately from Eq. ( |67] i that C k is the Fourier transform of (C" 4 (Ax)) = C(Ax), or one can plug 
Eq. ((68]) into Eq. ((65), to get (c£) = C k . 

E.2 Variance and Convergence 

The ergodic hypothesis is that an average with respect to a parameter such as position or time tends towards an 
ensemble average. In this case, 

w (Of) = C k , (69) 
and C A (Ax) « (C* A (Ax)) = C(Ax). 

when the surface area is very large. The ensemble average is a good substitute if the variance about the average 
vanishes as the substrate area A becomes large. It is found that in reciprocal space, 



Var(^)= -<C k 4 ) ^C*. (70) 
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Thus, the ergodic hypothesis does not hold for C£. In practice, is a speckled version of C k (Fig. [6j) However, if 
one smooths C£ by averaging over a small patch in reciprocal space of size fe smoo th = 1/A S , so that 

C£(A S ) = [§y 2 J d^e-^'-^C^ (71) 

then Var (C^ 1 (A s )) diminishes as 1/A. For sufficiently large A s , 

(C£(A S )) « C k , (72) 

and 

7r d/2 / \d 

Var(c£(A.))«^-^Cg. (73) 

Thus, the ergodic hypothesis (Eq. ( |69] i) only holds for a smoothed version of C^. 
In real space, 



Var(C^(Ax)) = ((^(Ax))] - (C A (Ax)) 
(2n) d 



A 



J d d k(e 2 * k * x Cl + CQ, (74) 



where the integral is bounded (finite) provided that either t > or the atomic scale cutoff 6q > 0. Thus, the ergodic 
hypothesis holds for the real space autocorrelation function. 

E.2.1 Eq. (JTOjl 

First, (C£C$) is calculated. 

A 



Assume that he distribution of /i k is gaussian. Also, assume that h(x) is real so that /ik^-k = |^k| 2 - Then, 

(/ikXAsO = CkxCk^^ki-k^^C^-kg)... 

... +C kl C k2 ^(k 1 + k 3 )^(k 2 + k 4 )... 

... +C kl C k3 <5 !i (k 1 -k 2 )^(k 3 -k 4 ). 

Thus, 



2 



2 



...+C k ^ d (k + k')] 2 + C k C k , [^(0)] 2 ). (75) 
= ^(<5 kk , + <5 k( _ k0 )+C k C k , ! (76) 
where Eq. d29| has been used liberally. Setting k = k', results in Eq. (|70|. 



E.2.2 Eq.(p3j 

Now consider C A smoothed over a length A s (Eq. (|7T|). The mean value is 

,2X^/2 



^k'e-lA^k'-k) 3 ^ 



2tt 
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For sufficiently small fc S mooth, (sufficiently large A s ), Eq. ( |72| results. 

The variance of (A s ) is now calculated. First, it is necessary to calculate / [C^(A S )] 2 



Using Eq. f75] l and Eq. |29) as needed, 



2n 



J d d k" e -s A ^ k "- k ) 2 (C&C&) 



2tt 



tfVrfV' e" 5 A * ( k '~ k ) 2 e" 3 A = ( k "- k ) 



••• x {c 2 [<J d (k' - k")] 2 + C 2 [6 d (k> + k")] 2 + C k C k , [^(0)] 2 } 
^ | d d k' {e- A ^ k '- k ) 2 C2, + e -5 A2 [( k '- k ) 2+ ( k ' +k ) 2 ]^,}... 



^1 

2tt 



2\ d/2 



d d k , e -lA^(k'-k) 2 c 



The first integral is bounded (finite) because C k is bounded. Let its finite value be denoted /. The second integral is 
simply (Cjf (A s )) .Thus, 



Var(C£(A.)) = 



Af J 



a finite value that decreases as A -1 as required for the ergodic hypothesis to hold. For sufficiently small fc smoo th (large 
A s ), I w (7r/Ag) d / 2 C^, and Eq. ( f73] > results. It should also be noted that the large A s required for this approximation 
also creates a more stringent requirement that A be large. 

E.2.3 Eq. (75) 

Now, consider the real space auto-correlation function. First, (C A (Ax)C A (Ax)) is needed. 

(C A (Ax)C A (Ax)) = J d d kd d k' e^+^'Ax (c£c£) 
Proceeding in a fashion similar to the previous section (making use of Eqs. fiT5\ and ([29]) as needed) , 



(C A (Ax)C A (Ax)) = /" d d kd d k! e '( k + k ')-Ax ^ C 2 ^ rf(k _ k ^ • 



Cl[6 d (k + k>)} 2 + C k C k , [5 d (0)Y 



A 



J d d k ( e 2 * kAx C 2 + Cl) 



^ d d ke lkAx C k j I J <fV e ik ' Ax C k 

(2n) d f „j, / 9il<-.Av ^.9 . ^9\ /^4/ • s\2 



.4 



J d d k ( e 2lkAx C 2 + Cl) + (C A (Ax)> 5 



Thus, Eq. ( f74| ) results. For the variance to be vanishing, the integral in Eq. fifty must be bounded (finite). If time, 
t > 0, the exponential in Eq. ( |77j ) guarantees that the integral is bounded. For time t = 0, the integral is only bounded 
if the atomic scale cutoff b > 0. 
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F Atomic Scale Cutoff 



Starting from Eq. d39k 



Ck= (2^F e2CTkMb ° fe2 - (77) 

The effect of the small scale cutoff is both small and short-lived, as it only works to suppress fluctuations with large 
wavenumbers. The most important fluctuations have wavenumbers between and 2k c . Thus, the typical size of the 
cutoff term is about b\k 2 c . If a typical dot size or spacing size 10 nm, and a typical atomic scale is 10 _1 nm, a typical 
value for this term is about 1CP 3 — 1CP 2 . To calculate the effect of the cutoff, it can absorbed into the time-dependent 
part with the substitution 

so that its effect lasts only as long as a perturbation with atomic scale curvature (k = 6 ). Thus, Eq. ( |40| i is a good 
approximation. 
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